\documentclass{article}
\usepackage{amsmath}
\usepackage{graphicx}

\begin{document}
\section{Shan-Chen model}
In Shan-Chen model, interaction forces between molecules are
accounted, and the interaction forces change the original ideal gas
equation of state into a non-ideal equation of state. With a
non-ideal equation of state, phase separation will take place when
the temprerature falls below the critical temperature of the fluid.
For multi component fulid simulation, intre-species interations are
also considered. \\

Compare with Naiver-Stokes momentum equations, the first term can be
included in the pressure term, and therefore results in a non-ideal
EOS as:

\begin{equation}
p=\rho c^2_s+\frac{1}{2}c^2_sG\Psi^2(\rho)
\end{equation}

The second term is related with interface stress corresponding to
the gradient free energy. It control the geometry of interface as
well as the surface tension. It is dependent on the mesh resolution
due to the existing of $\Delta x^2$.

In thermodynamic theory, the equilibrium pressures and densities
with particular temperature could be obtained by using Maxwell
construction (Rowlinson and Widom 1982). The Maxwell Construction
can be stated as
\begin{equation}
\int _{V_{m,l}}^{V_{m,v}} PdV_m = P{V_{m,v}-V_{m,l}}
\end{equation}

\begin{figure}[!hpg]
\begin{center}
\includegraphics[width=3in]{maxwellconstruction.eps}
\end{center}
\caption{The molar volumes of the liquid and vapor and vapor pressure can be determined by Maxwell Construction}\label{MxC}
\end{figure}



The area under the curve must equal the area of the rectangle
defined by vapor and molar volumes and the vapor pressure.Fig (\ref{MxC})\\

The Applying the mechanical equilibrium condition on the stress
tensor across the interface determine the interface structure.

\begin{equation}
\sigma= \int_{-\infty} ^{\infty} (P_{zz}-P_{xx})dz
\end{equation}

Assumption of interface lying in the x-y plane is made for above
equation.And $P_{zz}$ and $P_{xx}$ are the stresses normal and
parallel to the interface respectively. The value of surface tension
in the interaction potential  model have been reported by
Shan(2008)\cite{shan2008}:

\begin{equation}
\sigma=-\frac{eGc_s^4 e_4 \Delta x^2}{2} \int_{-\infty} ^{\infty}
(\frac{d\Psi}{dz})^2 dz
\end{equation}

From above analysis, we can find that, the surface tension is
related with the chosen of potential function which is equivalence
to the EOS, and as well as the lattice geometry. In order to make
the surface tension independent of  of EOS

By introducing a new term to instead of the second term in equation
(REFNUMBER) It's possible to make the surface tension independent
from the selection of EOS. Change the weight of the term $\psi
\nabla \nabla ^2 \psi$ with a coefficient $\kappa$

\begin{equation}
F^s = \kappa \psi \nabla \nabla ^2 \psi
\end{equation}

The discretization form of above function could be obtained by using
the pair-wise interaction between lattice site:

\begin{equation}
F^s (x) = \kappa \psi(x) \sum _{y} \lambda_{xy} \psi(y) |x-y| /
\Delta x^3
\label{addShanChenterm}
\end{equation}

A larger set of lattice sites in addition to the nearest neighbors
are need to obtain higher order difference operator for the
potential function $\psi$. The coefficient $\lambda_{xy}$ can be
determined by Taylor expansion of function
(\ref{addShanChenterm}).As a result, the total interaction force can
be written as:

\begin{equation}
F=F^0 + F^s \approx -\nabla (\frac{1}{2}c_s^2 G \Psi ^2) + \kappa
\Psi \nabla \nabla ^2 \Psi
\end{equation}




\section{MRT}
LBGK use single relaxation time and a Taylor series expansion of the
Maxwell-Boltzmann equilibrium distribution function for the
equilibria. There are some limitations for viscosity, low viscosity
will cause numerical instability. There is no conclusion for the
reason of this numerical instability, some have attributed them to
the nonexistence of an H theorem. (B.Boghosianl 2001)\cite{boghosianl2001}.
There exist two school of thought on how to overcome numerical
instabilities(McCracken 2005)\cite{mccracken2005} The first suggests using nonpolynomial
equilibrium distribution functions, however, doing this will
significantly increase the computing requirements and will not
always lead to correct Naiver Stokes equations. The others propose
the use of multiple relaxation times in order to reduce the
numerical instabilities. \\

MRT allows independent adjustment of bulk and shear viscosities.


\section{nonideal gases}
The continuous Boltzmann equation has mainly been used to solve
supersonic flows (K,Xu 1994)\cite{xu1994} This is due to the extreme complexity
of the collision kernel when dealing with dense, interacting
particles and the high requirement for computer resources. The
original Shan-Chen potential model only include attraction force
between particles, however, when the gas is compressed, the
repulsive force will dominate the gas model. Enskog theory take into
account the effects of molecular volume on molecular transport
properties.However, the long-range molecular interaction among
molecules is ceglected in Enskog's theory, we use mean-field theory
to describe the long-range molecular interactions. (He, Shan and
Doolen 1998) proposed a kinetic model that combines Enskog's theory
for dense fluids and mean-field theory for long-range molecular
interaction.

Continuum Boltzmann equation with BGK collision model is used to
built the kinetic model:

\begin{equation}
\frac{\partial f}{\partial t}+ \xi \cdot \nabla f+ F \cdot
\nabla_{\xi}f = -\frac{f-f^{eq}}{\lambda}
\end{equation}

Where $f(x,\xi,t)$ is the distribution function in the phase space
$(x, \xi)$, $\xi$ is the microscopic velocity and $F$ is the an
external body force. $f^{eq}$ is the Maxwell-Boltzmann distribution
function.

The macroscopic quantities are written as:

\begin{eqnarray}
&\rho=\int f d \xi& \\
&\rho \mathbf{u}= \int \xi fd\xi&\\
&\frac{D}{2} \rho RT= \int (\xi- \mathbf{u})^2 f d\xi&
\end{eqnarray}

The term $\nabla _{\xi} f$ cannot be computed directly since the
distribution function is dependent on the microscopic velocity which
is unknown. An approximation of using equilibrium distribution
function is made to compute this term:

\begin{equation}
\nabla _{\xi}f \approx \nabla_{\xi} f^{eq} = -\frac{\xi -
\mathbf{u}}{RT}f_{eq}
\end{equation}

Consequently, simplified Boltzmann equation is obtained as:

\begin{equation}
\frac{\partial f}{\partial t}+ \xi \cdot \nabla f=
-\frac{f-f^{eq}}{\lambda}+\frac{F \cdot (\xi -u)}{RT}f^{eq}
\label{equationenskog}
\end{equation}

There are two important factors for nonideal gas: the intermolecular
attraction and the exclusion volume of molecules. In He Shan and
Doolen model, the intermolecular attraction is treated using
mean-field approximation. Averaged force potential is used to
describe the intermolecular attraction force:

\begin{equation}
V(r_1)=\int _{r_{12}> \sigma} u_{attr}(r_{12})\rho(r_2)d r_2
\end{equation}

Where $r_{12}= |r_1-r_2|$ and $u_{attr}(r_{12})$ is the attractive
component of the intermolecular pairwise potential. $\sigma$ is the
diameter of the molecules. Expanding $\rho(r_2)$ about $r_1$ and
with the assumption that $\nabla \rho$ is small, approximation of
$V$ is achieved with first two leading term:

\begin{equation}
V=-2a\rho -\kappa \nabla ^2 \rho
\end{equation}

Where $a$ and $\kappa $ is constant given by:

\begin{eqnarray}
&a=-\frac{1}{2}\int _{r>\sigma} u_{attr} (r) d r&\\
&\kappa =\frac{1}{6} \int_{r>\sigma} r^2 u_{attr} (r) d r&
\end{eqnarray}

The effect of the exclusion volume of the molecules coule be
approximated by(S.Chapman and T.G.Cowling, the mathematical theory
of non-uniform gases 1970):

\begin{equation}
-f^{eq}bp\chi (\xi-u) \cdot \nabla ln(\rho^2 \chi)
\label{exclusionv}
\end{equation}

where $b=2\pi \sigma^3/3m$ $m$ is the mass of a single particle.
$\chi$ is the increase in collision probability due to the increase
in fluid density:

\begin{equation}
\chi = 1+\frac{5}{8}b\rho + 0.2869(b\rho)^2 +0.1103 (b\rho)^3+
\ldots
\end{equation}

Combine the equation (\ref{exclusionv}) with the intermolecular
attraction term, a equivalent force term is obtained:

\begin{equation}
F=-\nabla V-b\rho RT\chi \nabla ln(\rho^2 \chi)
\end{equation}

Chapman-Enskog analysis yield the corresponding macroscopic momentum
equation as:

\begin{equation}
\frac{\partial (\rho u)}{\partial t}+ \nabla \cdot (\rho uu)=-
\nabla \cdot P' +\nabla \cdot \tau + \rho a
\end{equation}

Where $\tau$ is the viscous stress $\tau=\mu_E S$ viscosity $\mu_E$
is defined as:

\begin{equation}
\mu_E = \rho c_s^2 (\frac{\tau}{\chi}-\frac{1}{2})\delta_t
\end{equation}

$P'$ is the thermodynamic pressure tensor with the definition as:
\begin{equation}
P'_{\alpha \beta} = p\delta_{\alpha \beta}+\kappa \frac{\partial
\rho}{\partial x_{\alpha}}\frac{\partial \rho}{\partial x_{\beta}}
\end{equation}
\begin{equation}
p=\rho RT(1+b\rho \chi)-a \rho^2 -\kappa \rho \nabla^2 \rho
-\frac{\kappa}{2} |\nabla \rho|^2
\end{equation}

The term $\nabla \cdot P'$ can be written as:
\begin{equation}
-\nabla \cdot P'=-\nabla p_0 +\kappa \rho \nabla \nabla^2 \rho
\end{equation}

The first term $\nabla p_0$ on the right hand side is the gradient
of equation of state:

\begin{equation}
p_0=\rho RT(1+b\rho \chi)-a \rho^2
\end{equation}

The second term $\kappa \rho \nabla \nabla^2 \rho$ represent the
surface tension. Different equation of state can be employed by
modifying the term value $p_0$.The force term can be written as:

\begin{equation}
F_t=\rho a_t = -\Delta \psi +F_s+\rho a \label{forcetermenskog}
\end{equation}

$\psi$ is a function related with equation of state:
\begin{equation}
\psi=p_0-p_{id}=p_0-\rho RT
\end{equation}


To solve equation \ref{equationenskog} numerically, we first
discrete the equation in time:
\begin{equation}
f(x+\xi \delta t,\xi, t+\delta t)-f(x,\xi,t)=-\int _t^{t+\delta t}
\Lambda (f-f^{eq}) dt+\int _t^{t+\delta t} \frac{F \cdot
(\xi-u)}{RT}f^{eq} dt
\end{equation}

The first integration is assumed to be constant over one time step,
and a trapezoidal rule is used to compute the integrand of the
second term

\begin{equation}
f(x+\xi \delta t,\xi, t+\delta t)-f(x,\xi,t)=-\Lambda (f-f^{eq})
\Delta t+ \frac{1}{2}(\frac{F \cdot (\xi -u)}{RT}f^{eq}|_{t}+\frac{F
\cdot (\xi -u)}{RT}f^{eq}|_{t+\Delta t})
\end{equation}




The LBE in the general form is:

\begin{equation}
f_i(x_i+c_i\Delta t, t+\Delta t)-f_i(x_i,t)= -\sum_j
\Lambda_{ij}(f_j-f_j^{eq})|_{(x,t)}+\frac{1}{2}[S_i
(x,t)+S_i(x+c_i\Delta t, t+\Delta t)]\label{MRTenskog}
\end{equation}

where $c_i$ is the discrete velocity, $S_i$ represent the force term
in Eq. (\ref{forcetermenskog}) The equation (\ref{MRTenskog}) has a
implicit term which usually requires inverse calculation. To
maintain the explicit scheme, we need eliminate this implicity, an
new variable is introduced to do so

\begin{equation}
\bar{f}_i =f_i-\frac{1}{2}S_i \Delta t \label{newtranformation}
\end{equation}

With the transformation in equation (\ref{newtranformation}),
equation (\ref{MRTenskog}) can be converted into an explicit scheme
as:

\begin{equation}
\bar{f}_i(x+c_i\Delta t, t+\Delta t)-\bar{f}_i(x,t)= -\sum_j
\Lambda_{ij} (\bar{f}_j-f^{eq}_j)|_{(x,t)}+\sum_j
(I_{ij}-\frac{1}{2}\Lambda_{ij})S_j|_{(x,t)}\Delta t
\end{equation}

Where $I$ is an identity matrix. The kinematic viscosity $\nu$ is
defined as
\begin{equation}
\nu=(\frac{1}{s_8}-\frac{1}{2})c_s ^2 \Delta t
\end{equation}

And the bulk viscosity is related with $s_2$:
\begin{equation}
\zeta =(\frac{1}{s_2}-\frac{1}{2})c_s^2 \Delta t
\end{equation}
 The macroscopic quantities can be calculated from the summations:
\begin{eqnarray}
&\rho=\sum_{\alpha} f_{\alpha}=\sum_{\alpha} \bar{f}_{\alpha}&\\
&\rho u=\sum_{\alpha} f_{\alpha}e_{\alpha}=\sum_{\alpha}
\bar{f}_{\alpha}e_{\alpha}+\frac{1}{2}F\delta t&
\end{eqnarray}





\section{Literature review}
Due to the interaction are limited with attraction force only, when
the gas is compressed to near its 'hard sphere' volume, the
repulsive force will dominate the van der Waals gas model.
Thermodynamic theorem is satisfied only when the potential function
is in the form of exponential function:

\begin{equation}
\Psi (\rho)=\psi e^{-\rho_0/ \rho}
\end{equation}

Therefore the original Shan-Chen model has very limited ability to
simulate any desired EOS. Incorporating with Enskog dens gas theorem
can resolve EOS limitation. Several LBM incorporate with Enskog dens
gas theorem has been developed. (He and Doolen 2002, Martys 2001,
Luo 2000, Luo 1998)\cite{he2002} \cite{martys2001} \cite{luo2000} \cite{luo1998}. 

A discrete model based on the Boltzmann equation with a body force was developed by He,Shan and Doolen(1998)\cite{he1998}. The Boltzmann equation is discretized in a way that preserves the derivation of the hydrodynamic equation from the Boltzmann equation. A mean-field approximation is used in their model for the interparticle interaction.They reported that, in this type of model, anisotropy might be inevitable if a nonlocal interaction is not included in momentum space. 

LBM multiphase models based on Enskog theorem have been extended with MRT algorithm in 2D and 3D(McCracken,2005, Premnath, 2008)\cite{mccracken2005} \cite{premnath2008}. Improvement in stability and capability for simulating low viscosity fluid was reported to be achieved.

Comparison of different multiphase LBM model in thermodynamic aspect has been given by He and Doolen (2002)\cite{he2002}. Detailed analysis of thermodynamic foundations of kinetic theory was given. Demonstration shows that kinetic model is consistent with thermodynamic theory and can lead to correct not only mass and momentum equations but also energy transfer equation.

The spurious problem was first studied by Shan (2006) \cite{shan2006}. Analysis show that the spurious current near interface is due to the insufficient isotropy of the discrete gradient operator. Finite difference gradient operators with higher order of isotropy was proposed and the spurious current was found decreasing significantly.  Sbragaglia (2007)\cite{sbragaglia2007} analysis the flows near interface in details and extended pseudopotential method is developed with using multirange pseudopotential. The extended pseudopotential model also permits to tune the equation of state and surface tension independently of each other.



\bibliographystyle{ieeetr}
\bibliography{ref}

\end{document}
